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The presented path to PDFs extraction, via global analysis of data, is based on preassumed form 
of parametric PDFs. In this paper, PDFs are obtained at a chosen Q 2 , and at Xi,i — 1, n of the 
analyzed data points by solving the numerically singular system of decomposition equations of the 
data points each at a given {xi,Qij 2 ) of the data. In each equation, the corresponding data point 
is decomposed into its pdf components, evolved from the set of unknowns at (xi, Q 2 ), i = 1, n to 



the required (xi, Qi 3 
evolution equations. 



,i = l, 



,n,j = 1, 



via our linear algebraic (matrix) solution to DGLAP 
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INTRODUCTION 

Given conditions on parton distribution functions 
(PDFs), coming from existing measurements in the accel- 
erators, e.g. from a measurement of the structure func- 
tion F 2 in DIS, through a fixed order (LO, NLO, ... ) 
formula from PQCD for F%, or a given cross section for pp 
scattering in a similar manner, can be considered as the 
initial conditions for the PDFs in the integro-differcntial 
DGLAP evolution equations of the same order in pertur- 
bation; based on which solving the DGLAP equations 
generally allows PQCD prediction for processes of in- 
terest in accelerators. Thus, the method of solution of 
DGLAP proposed in the year 2000 Q, and developed here 
can be of serious interest. For simplicity the method is 
illustrated at LO. 

DGLAP evolution equation is of the simple one dimen- 
sional form for the non-singlet quark distribution /ns'- 



Where 



df N s(x,t) f dy x 

=/ —Pqq(-)fNs(y,t) 

ot J x y y 



33-2n fr / ,Q 2 X 



(1) 



(2) 



In the more traditional form of DGLAP the L.H.S. 
derivative is with respect to Ln(Q 2 ). 



2n.dt 

dLn(Q 2 ) 



a s {Q 2 ) 
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(33 - 2n f )Ln{%) 



(3) 



Here, the strong coupling constant a s (Q 2 ) is rendered in 
LO in the historical approach^, using a constant of in- 
tegration A 2 = A 



QCD 



being in the range 200 — iOOMev, 
signaling the QCD confinement, where the coupling be- 
comes very strong; e.g. as hinted in LO, ©. We will use 
the LO formula ([3]) for the strong coupling constant con- 
sistently throughout the calculations later on, avoiding 
the ambiguities and resulting mistakes that may arise if 
definition is varied. The present day approach bypasses 



the ambiguities by directly using experimental determi- 
nation of a s {M z ) with no recourse to Aq C£) [^, 

DGLAP evolution equation has two dimensional form 
for the singlet quark distribution E, coupled to gluons: 
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dt 



T,(x,t) 
g(x,t) 



'dyfPqqi^) 2n/P TO (|: 
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At LO, where equality of splitting functions (neglecting 
quark masses), leaves only four different ones, the evo- 
lution equation can be written in a simple (2n + 1) di- 
mensional form for the rif quark and ant-iquark flavors 
coupled through the gluon distribution: 
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ot 



(l2n s (X, t) 



dy 



(Pqq(^) 






P, 
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P M (|)P, 9 (|) 
Pgq(y) Pgg(y)/ 



(5) 

( Qi{v,t) \ 

Q2n f (V,t) 

V g(y,t) ) 



N-DIMENSIONAL DISCRETIZATION OF 
X-SPACE OF THE DATA POINTS AND PDFS. 

Analytical solutions to DGLAP exploits nearness of 
the discrete data points of a given set, e.g. of the struc- 
ture functions, to the Bjorken Xi, 



Xi = (xiY, 

< ... < Xi + i < Xi < ... < x\ < l,i = 1, n. 



(6) 
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Numerical solutions can simply use the exact Bjorkcn 
Xi,i = l,...,n, of the data set. In either case, a discrete 
set of n basis vectors for the PDFs is defined on each 
Xj. Thus, the splitting functions operating on this basis 
(space) can be calculated asnxn matrices. 

The discrete basis for the PDFs are used in the con- 
volution integrals on the right side of the DGLAP equa- 
tion to produce the matrix form of the kernel K, which 
in DIS data analysis, is either the splitting functions 
for DGLAP, or the coefficient functions of the hadronic 
structure function. The convolution integrals are of the 
generic form: 



The convolution ([7]) becomes: 



P®f 



1 ^L K {*)f{y,Q 2 ) 

v y 



1 ^K{ y )f{*,Q% 

y y 



(7) 



in which, for the simplest parton distribution function, 
e.g. non-singlet (NS) f(x,Q 2 ) = f NS (x,Q 2 ), the 
discrete basis is an n— tuple, whose ith element is 
f(xi,Q 2 ). As the convolution integral is made a discrete 
sum, Pij,i, j = q,g or Ck,k = q,g within the kernel 
K become lower triangular (l.t.), due to the limits of 
integration of the resulting sum; and banded (bdd), 
due to divisions of Xi in fl}. The resulting algebraic 
simplicity of the l.t. bdd matrices is an essential element 
for having analytical solution of DGLAP as the basis of 
the analysis. Evidently, numerical solution working with 
exact Xi,i = 1, ...,n of the data neither works with, nor 
requires handedness. 
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Afl + y 2 



i-y 



(10) 



where the extra y in the left integrand is put in to work 
with the more useful form, the momentum distribution. 
Note that the DGLAP equation can be written for mo- 
mentum distribution, i.e, in f(x) —> xf(x), by putting 

an extra x on both sides, and a factor of v.- within the 

. v . 

integral. Different techniques may be used for doing the 

integration in DGLAP. Here, we use only integration by 

parts: 

u = /(-) ^>du = — — dy, dv = -(- )dy =>• v = / dv; 

y dy 3 1 - y J 

So, 

h = f(x)v(l) f(l)v(x) J' v(y)^dy =f(x)v(l)+ 

1 ,x df(x) 
y dy 

l 2 = -f(x)[v(l)-v(0)}. 



(11) 



Illustration of the method by an example 



fv{-)^P-dy. 
y dy 



To calculate the matrix form of a splitting function, 
e.g. P qq in ((T|), a combination of discrete finite differ- 
ence approximation of the unknown PDFs and continu- 
ous integration of the known splitting functions is mixed 
in our method of evaluation of the convolution, as fol- 
lows. Then, (Pqq)ik is extracted from the coefficients of 
fk = f(xk,t). At LO, P qq is given in the following equiv- 
alent forms |5I1 : 



-1 
M 3 



1 



y 



i-y 



l+V 

(1-2/)- 



+ b(l-y) 



(8) 



Given the defnition of "+" rcgularization, aimed at re- 
moval of the infinity at y = 1 in the kernel P qq in the 
context of the convolution integral^: 



dxf(x) 



l-x 



(9) 



l-x 



g(x) - /(l) J dxj 



Here and later, when convenient, the variable change 
^ — > y is used. Bringing in a finite difference approx- 
imation of the differentials: 



I = v(0)f( Xi ) + J2 { v(xi/y) 

k=l Z 



fk — fk- 
Xk - X k -1 



-dy 



w(0)/(xi)+ / v(x t /y) fi fi ' dy 



i-i -"r 1 

/ v(xi/y) 



fk - f, 



fc-1 



Xk - Xk-i 



dy 



v(0)f(x t )-\ / v(x l /y)dy- 

***% 



fi- 



- — / u{xi/y)dy + ^2 / vfe/y) 

i—1 J k—1 



(12) 



fk — fk- 
Xk - x k -i 



-dy. 
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Now, a bdd l.t. P qq can be read (extracted) from the 
coefficients of ff. . 



P 



(P qq )u = «(0) + 
1 

(Pqq)ik 



1 



v(xi/y)dy, 



x k - x k -i 



v( Xl /y)dy- (13) 



v{x l /y)dy. 



Xk + l 

The other three kernels are essentially obtained in a sim- 
ilar manner. In addition, construction of the matrix co- 
efficient functions in MS which are taken fromQ, follows 
closely the path of construction of splitting functions, re- 
sulting in commuting, bdd, l.t. or non-commuting, non- 
bdd, l.t. matrices, respectively, in analytical or more 
general numerical approaches. 



LINEAR ALGEBRAIC (MATRIX) 
FOUNDATIONS OF THE NUMERICAL 
SOLUTION OF DGLAP IN X-SPACE 

In this section, the numerical solution of the DGLAP 
equations (JTJ, Q and (|SJ) is explained in detail. The 
DGLAP equations are of the generic form: 



df K (x,t) 
dt 



1 <^K{ X -)f K {y,t), 
V V 



(14) 



Where K = Pns,Ps,Pf is one dimensional in non- 
singlet, two dimensional in singlet and (2rif + 1) dimen- 
sional in flavor kernels, respectively. 
Having the matrix forms of the above kernels from sec- 
tion (2), the matrix form of DGLAP equations reads as: 

^^ = E^/f(i)^T X (t)=K.f^(t); (15) 

Since the kernel K being independent of the variable t, 
leads the following solution to the matrix equation (fTSI ) 
for finite evolution from to to t: 



f K (t) = e^-'^.f^to) = Ef .{ K (t 



(16) 



Therefore, the finite evolution machinery for the DGLAP, 
is constructed with the final (t — to) dependent Ef ma- 
trix. 

To determine the matrix form of Ef = e K( -* ~^ in 
analytical framework of [l[ , as mentioned in the previous 
sections, the choice of x% = x\, i = 1, • • • , n is manda- 
tory to have commuting bdd splitting functions matrices, 
but the experimental data are not always match these x 
points, so interpolation techniques should be used to shift 



the data points to the desired places. 
On the other hand, given the numeric form of the LO 
kernels PnSi Ps an d Pp to a software such as Math- 
ematica (using the MatrixExp command), immediately 
gives the matrix form of the exponential in the solution 
(|16l) . Using numerical algorithms to approximate the in- 
finite series of the exponential has the advantage of being 
needless of the commuting bdd matrices, so, ones could 
chose the experimental x to discrete the x space. 



COMPARISON OF THE EVOLUTION MACHINE 
RESULTS WITH MSTW AND CTEQ 

The machinery thus far has brought us to the point 
that we may take a discrete set of input PDFs at x-points 
patterned according to ((BJ) at to from some present stan- 
dard LO set, such as MSTW2008, evolve it ourselves to 
t, and then compare our results with their PDFs at t. 
This is actually done via an n = 100 (O input taken 
at Qq = lOGev 2 from MSTW (green curves). Evolution 
needs the value of a s which, in the line of the discussion 
of equation ([3]), was handled using Aqcd consistently. 
The excellent quality of match of our solutions outputs 
at Q 2 = lOOGew 2 (red curves) and MSTW PDFs (blue 
curves) can be seen in Fig.(jT]-I) for non-singlet - singlet, 
and in Fig.([T]-II) for the flavors. 

Definitions of 5 flavor non-singlets and the singlet are: 
q? S = (u + u)-(d + d), 



NS = (u + u) + {d + d)-2(s + s), 



q^ s = {u + u) + {d + d) + (s + S) - 3(c + c), (17) 



q^ s = (u + u) + (d + d) + (s + s) + ( c + c) - 4(6 + b), 
£ = (u + u) + {d + d) + (s + s) + (c + c) + (b + b). 

A surprise is that, we assume complete SUf (5) (flavor) 
symmetry with our up to , or flavor DGLAP Kernel 
([5])- Then, an n = 100 point fit has such a good match 
with MSTW2008 that says they should be having such an 
SUf (5) symmetry assumption as well. The only source 
of asymmetry in both our procedures seems to be the 
non- symmetric PDFs inputs to the totally symmetric 
evolution process, resulting in the asymmetric outputs. 
A similar comparison is shown in figurc([T]-III) with the 
older LO CTEQ5L set. It is interesting to note that the 
excellent match, points to their SUf (5) symmetry as- 
sumption as well. 

There is only a final step needed to have the detailed 
match observed in the figures. While solving DGLAP, 
flavor number changes should be accommodated every 
time a changing Q 2 crosses a massive quarks m 2 . Ac- 
cordingly, a simple method of continuous evolution across 
stages with breaks at Q 2 = m 2 was conceived, which may 
be called zero-mass scheme. A more involved zero-mass 
scheme is found in MSTW But there is no apparent 
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difference for our LO comparison! In both, as mentioned 
there, the name zero-mass is misleading, as mass appears 
in the boundary conditions. Continuity of PDFs is kept 
via the boundary conditions of the DGLAP at breaks at 
Q 2 = m 2 ,. Specifically, a break takes place if evolution of 
the PDFs, during increase of t(Q 2 ) ([3]), from Uf flavors 
at to, crosses a t q = t(m q ), of a massive quark h, to reach 
t with rif+1 flavors. Then, solving DGLAP for finite evo- 
lutions [to,t] is divided into two stages of [to, t q ] &nd[t q , t\. 
Furthermore, there are two sets of inputs, the first set has 
n; flavors at to, the second set has (n/ + 1) flavors at t q . 
In addition, the second set of inputs is the first set of 
outputs with an extra input PDF for the heavy quark, 
h(m 2 ) = 0. Specifically, in our figures, h= b-quark with 
MSTW m b = A.75Gev, and CTEQ m b = 4.5Gev. 



EXTRACTION OF PDFS FROM F 2 p - d DATA 

Beginning with LO structure function F2 symbolically 
as: 

F2 = X{1 + 2^° q) ® P2+ X ^n° 9 ® E e /)5- (18) 

with the quark structure function in term of rif = 5 fla- 
vors NS-S, with usual definitions ([P?])^ , 

F 2 = ^( 22S + 3? (24) _ 5 <?(15) + 5q {8 ) + 15g (3) ), 

(19) 

F2 1 = 2 2 2 = ^(22S + 3<z (24) - 5q m + 5g (8) ). 

Now, we are ready to actualize the program explained 
at the onset. Each F^' d data point at a given (xj,Qf-) 
of an experiment, e.g. BCDMS0, is decomposed into 
its partonic components via (TTg)) and (TT5]), evolved by 
our solutions of DGLAP, from a set of unknown parton 
distributions at (xi,Q 2 ) to be solved for by the linear 
system of equations for all the data and the variable flavor 
number (VFN) constraints. Our data equations come 
from minimization of x 2 with respect to the unknowns, 
with the usually defined 

X 2 = £ J ( 20 ) 

i=l 1 

where F^ is the value of data point on the left side of (|T8|) 
with its quadratically calculated total error being cr^, and 
F^ is its LO decomposition corresponding to the right 
side of (fl8|) . For the range of Q 2 used, the VFN con- 
straint equations are set at the masses of c and b quarks, 
considered as boundary conditions (b.c.) between regions 
of evolution in Q 2 where tif respectively crosses 3- 4, and 
4-5. Constraints equate the singlet, S, and the (n 2 f — 1) 



-multiplct NS of the SU{nf) flavor symmetry above each 
b.c; thus, dropping the flavor number, nt, by 1 at the 
b.c. The constraint equations, together with consistent 
evolution, keep the continuity of the PDFs while the fla- 
vor numbers changes at the b.c. 



RESULTS 

As mentioned in the previous section, each data point 
serves as a constraint equation on PDFs. A measurement 
of the structure function F2 in DIS, through a fixed or- 
der (LO, NLO,. . .) formula from PQCD for F2, or a given 
cross section for pp scattering in a similar manner, can 
be considered as such a constraint equation. The integro- 
differential DGLAP evolution equations of the same or- 
der in perturbation may then be solved and used to tic 
all the data points in an over determined and numerically 
singular system of linear equations to be solved via sin- 
gular value decomposition (SVD)[iJ together with other 
constraint equations on the PDFs coming from variable 
flavor numbers (VFN) pertinent to the Q 2 range of the 
data. Here, the presented method of solution of DGLAP, 
becomes of serious interest. The presentation here is done 
at LO, for the simplest example of DIS data, F%' d from 
BCDMS collaboration^ ; however, it is treated in gen- 
eral terms for other applications and applicable to be 
extended through PDF Global Data Analysis as well as 
higher orders of perturbation. 

Solutions for the PDFs are presented in figure©. The 
error analysis comes from error propagation equation (lo|. 
applied on the output of the solved system of equations. 
For simplicity, we assumed no correlation among experi- 
mental errors. Total number of used data points, for pro- 
ton and deuteron respectively, are n p = 112 and rid = 96. 
Number of data equations (same as the number of un- 
knowns) obtained from \ 2 are 6(PDFs) x ll(x-points), 
and those obtained from b.c. are 2(b.c.) x ll(x-points) 
and x 2 /N pts . = 0.93. Noting that MSTW2008LO PDFs 
are coming from global fits with the total number of 2615 
data points from 36 data sets. 



CONCLUSION 

We have developed in detail, analytical linear alge- 
braic solution to DGLAP, with its numerical presenta- 
tion, employed it in doing non-parametric data analysis; 
i.e. obtaining PDFs directly from data points without a 
pre-assumed parametric form. Our results for finite LO 
evolution match those of the simply available standard 
LO sets: MSTW2008 and CTEQ5L. In addition, the ob- 
tained parton distributions agree with the sophisticated 




FIG. 1. (a-I, II, III) LO Parton distributions extracted from MSTW2008 Ns-S, Flv. And CTEQ5L Ns-S respectively at 
Q 2 = WGev 2 . (b-I, II, III) Parton distributions at Q 2 = lOOGeu 2 (blue curves) in comparison with its corresponding evolved 
distributions from Q 2 , = WGev 2 to Q 2 = lQQGev 2 using the mentioned technique (red dashed curves). 




FIG. 2. Reasonable matches among extracted PDFs (red dots) from BCDMS F^ d DIS data in the range of x = 0.07 to x = 0.75 
in comparison with MSTW2008LO global fit distributions (blue curves) at fixed Q 2 = 37.5Get; 2 . 
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MSTW and CTEQ global data analysis. 
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